L  POSTGRADUATE  SCHOOL 

Monterey,  California 


FINITE 

METHODS 

FOR 

A 

NONLINEAR  ALLOCATION 

PROBLEM 

ALAN 

R.  WASHBURN 

APRIL  1989 

Approved   for  public   release;    distribution   is   unlimited. 


PedDocs  ared  for: 

D  208.14/2 

NPS- 55-89-003   ,f  of  Naval   Operations 
lington,   D.C. 


SOBOOI. 


NAVAL  POSTGRADUATE  SCHOOL, 
MONTEREY,  CALIFORNIA 


Rear  Admiral  R.  C.  Austin  Harrison  Shull 

Superintendent  Provost 


The  work  reported  herein  was  supported  in  part  with  funds  provided  from  the  Chief 
of  Naval  Operations. 

Reproduction  of  all  or  part  of  this  report  is  authorized. 

This  report  was  prepared  by: 


iCURITY  CLASSIFICATION  O^  THiS   PAGE 


DUDLEY  KNOX  LIBRARY 

HAVAL  POaTOFlADUMI  b  bUHUOL 


REPORT  DOCUMENTATION  PAGftONTEREY  CA  93943-5101 


a    REPORT  SECURITY  CLASSlF,CATiON 
UNCLASSIFIED 


b    RESTRICTIVE   MARKINGS 


a.  SECURITY  CLASSIFICATION  AUTHOc 


b    DECLASSIFICATION    DOWNGRADING   SCHEDULE 


Distribution  ava   .-        t  of  repor' 
Approved    for   public    release;    distribution 
is   unlimited. 


PERFORMING  ORGANIZATION  REPORT  NUMBER(S) 

NPS55-89-03 


5    MONITORING  ORGANIZATION   REPORT   NuMBERi'S, 


a    NAME  OF  PERFORMING  ORGANIZATION 

Naval   Postgraduate    School 


6d    OFFICE  SYMBOL 
(If  applicable) 

Code   tf 


7a    NAME  OF   MONITORING  ORGANlZA'.ON 


c.  ADDRESS  {City,  State,  and  ZIP  Code) 

Monterey,    CA        93943-5000 


7o    ADDRESS  (City.  State,  and  ZIP  Code) 


a.  NAME  OF  FUNDING/ SPONSORING 
ORGANIZATION 

Naval   Postgraduate  School 


Jb    OFFiCE  SYMBOL 
(If  applicable) 


9    PROCUREMENT  INSTRUMENT   IDENTIFICATION   NUMBE: 

0  &  MN,   Direct  Funding 


c.  ADDRESS  (City,  State,  and  ZIP  Code) 

Monterey,   CA.     93943 


10    SOURCE  OF  FUNDING   NUMBERS 


PROGRAM 
ELEMENT  NO 


PROJECT 
NO 


■ 
NO 


■j  N I T 
ACCESSION   NO 


1.  TITLE  (Include  Security  Classification) 

Finite  Methods  for  a  Nonlinear  Allocation  Problem 


2    PERSONAL  AUTHOR(S) 

Alan   R.    Washburn 


3a.  TYPE  OF  REPORT 
technical 


13b    T:ME  COVERED 

from    Oct    87    to   Sep 


4    DATE  OF  REPORT    (Year,  Month,  Day) 

April   1989 


5    PAGE  CC 
20 


6.  SUPPLEMENTARY  NOTATION 


COSATi  CODES 


FIELD 


GROUP 


SUB-GROUF 


18    SUBJECT  TERMS  (Continue  on  reverse  if  necessary  and  identify  by  block  number) 
allocation,    nonlinear 


9    ABSTRACT  (Continue  on  reverse  if  necessary  and  identify   by  block  number) 

We  consider  the  problem  of  allocating  m  resources  to  n  tasks.   The  potential  effective- 
ness, y.  for  the  j    task  is  a  linear  function  of  the  allocations  to  that  task,  and  the 
objective  function  is  a  convex  function  f  (y  ,  . . . ,  y  )  to  be  minimized.   Two  finite 
algorithms  for  obtaining  an  optimal  solution  are  proposed  and  evaluated. 


20    DISTRIBUTION /AVAILABILITY  OF  ABSTRACT 

3  UNCLASSIFIED/UNLIMITED       □  SAME  AS  RPT  □  QTIC  USERS 


21     ABSTRACT  SECURITY   CLASSIFICATION 

Unclassified 


!2a    NAME  Oc  RESPONSIBLE  INDIVIDUAL 

A.    R.    Washburn 


22b  TELEPHONE  (Include  Area  Code) 

(408)    646-3127 


22c    OFFICE   5*   . 

55Ws 


)DFORM  1473,  84  mar 


APR  ed.tion  may  be  used  until  exhausted 
ther  editions  are  obsolete 


SECoR'TY  CLASS  FlCA'iON  C  AGE 

XI  U.S.  Govrnment   Printing  Office     1966—606-243 


Finite  Methods  for  a  Nonlinear  Allocation  Problem 

Alan  Washburn 


1.    Introduction 

We  consider  here  finite  algorithms  for  solving  the  nonlinear  program  PI : 


minimize  f(yi,...,yn) 

m 

subject  to  £  E(i,j)xij=yj;      l<j<n, 


i=l 

n 


£  xij=bi;         l<i<m, 

xij>0;  l<i<m,  l<j<n. 

It  is  assumed  in  PI  that  bj>0  for  all  i  and  that  E(i,j)>0  for  all  (i,j),  with  at  least  one  positive 
E(i,j)  in  each  row  i.  It  is  further  assumed  that  f(y)  is  continuous  and  decreasing  in  each  of 
its  arguments  on  some  convex  set  S  that  includes  the  nonnegative  orthant  of  Rn.  y 
represents  the  vector  (y  i,. .  .,yn)>  and  similarly  x  will  represent  the  collection  of  all  mn  of 
the  X]j,  etc.  PI  is  a  special  case  of  the  nonlinear  network  problem,  a  problem  that  is  already 
known  to  possess  sufficient  special  structure  to  make  the  development  of  special  purpose 
software  attractive(Ahlfeld,  Dembo,  Mulvey,  and  Zenios(1987)).  Our  intention  is  to 
specialize  PI  even  further  by  making  restrictive  assumptions  about  f(y).  The  components 
of  y  will  be  called  "potentials",  one  for  each  column  of  the  matrix  (E(i,j)).  Potentials  are 
not  logically  necessary  in  defining  PI,  since  the  expressions  defining  y  could  simply  be 
substituted  into  f(y),  but  the  generous  notation  will  prove  useful  in  the  sequel.  The 
components  of  x  will  be  called  "allocations." 

PI  could  be  interpreted  as  a  Search  Theory  problem  by  assuming  that  search  is  to  be 
conducted  over  a  fixed  time  period  by  m  distinct  types  of  searcher,  with  bt  units  of  type  i 
effort  available  over  the  period.  Assume  that  a  single  target  is  located  with  probability  pj  in 
one  of  n  regions,  that  region  j  has  area  Aj,  and  that  the  reward  for  finding  the  target  in 
region  j  is  Vj,  Assume  further  that  the  amount  of  area  swept  by  a  searcher  of  type  i  in 
region  j  is  Sjj.  If  xjj  is  the  allocation  of  search  effort  of  type  i  to  region  j,  then  the  total  area 

m 

swept  in  region  j  is  ^  XijSjj.  If  the  search  in  each  area  is  "random"  (Koopman(1956)), 

i=l 

then  the  number  of  times  a  target  in  region  j  is  detected  is  a  Poisson  random  variable  with 
mean   ^  xijsij  f  Aj  =yj  (so  let  E(i,j)=Sjj/Aj),  and  the  probability  of  no  detections  is 

n 

exp(-yj).  The  average  "regret"  in  terms  of  value  not  found  is  then  f(y)=2^  VjPjexp(-yj). 

j=i 
PI  is  then  the  problem  of  allocating  m  resources  to  n  regions  in  order  to  minimize  the 
average  undetected  value.  PI  can  be  thought  of  as  a  generalization  from  the  case  m=l, 
which  already  has  specially  tailored  solution  procedures  (Washburn(1981)). 

The  US  Air  Force's  Heavy  Attack  model  (Clasen,  Graves,  and  Lu(1974))  is  also 
closely  related  to  PI ,  with  resources  being  aircraft  sorties,  areas  being  classes  of  targets, 
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and  potentials  being  "average  number  of  targets  killed  were  it  not  for  the  fact  that  some  of 
the  potential  will  be  'wasted'  in  attacking  targets  already  dead".  In  general,  PI  is  applicable 
when  n  projects  must  be  simultaneously  undertaken  by  employing  m  resources,  y  is  a 
vector  measure  of  progress,  and  f(y)  is  some  global,  scalar  measure  of  work  not 
completed. 

2.    Basic   Feasible   Solutions 

PI  has  mn  potentially  positive  allocations,  but  optimal  solutions  invariably  have  most 
of  these  being  0.  The  reason  for  this  is  that  optimal  solutions  need  never  include  cycles  of 
positive  allocations. 

Definition:  If  a  feasible  allocation  for  PI  has  an  alternating  sequence 
Jl»il02»i2»--MiLJL+l  of  columns  and  rows  all  distinct  except  that  JL+i=jb  andif  xjj>0  for 
every  consecutive  pair  in  the  sequence  ((i,j)=(ik.jk)  or  (i,j)=(ik,  jk+l)  for  some  k<L),  then 
the  sequence  is  a  cycle. 

Definition:   A  feasible  allocation  for  PI  is  conservative  if  xjj=0  whenever  E(i,j)=0. 

Theorem  1:  Given  the  assumptions  of  section  1,  PI  has  an  optimal,  conservative 
solution  with  no  cycles. 

Proof:    PI  has  an  optimal  solution  because  the  objective  function  regarded  as  a  function 
of  x  is  continuous  on  a  compact  set  (Bazaraa  and  Shetty(1979)).  An  optimal  solution  x 
can  easily  be  converted  to  a  conservative  optimal  solution  by  shifting  any  offending 
allocation  xj ;  from  column  j  to  some  other  column  k  where  E(i,k)>0  (by  assumption  at 
least  one  such  column  k  always  exists  for  every  i).  Assume  therefore  that  x  is  optimal  and 
conservative,  but  that  a  cycle  exists.  In  that  case  construct  a  modified  solution  x'  by  adding 

8(i,j)  to  xjj  for  each  consecutive  pair  (i,j)  in  the  sequence,  where  5(i,j)  is  defined  as 
follows: 

8(ii,ji)=z  (an  arbitrary  real  number) 

5(ik0k+l)=-5(ik,jk);  l<k<L 

S(ik+lJk+l)=-  &(ik,Jk+l)E(ikOk+l)/E(ik+lJk+l);  13c£L-l. 

n 

The  numbers  5(i,j)  are  constructed  so  that  ^  x'ij=Di  for  all  i  as  before.  Also 

j=i 

m 

2^  x'jjE(i,j)=yj  as  before,  except  that  the  latter  sum  for  j=ji  is  increased  by 

i=l 

5(i i ,j  i )E(i i,j  i)-5(iL,j i)E(iL,j  i).  This  difference  is  proportional  to  z.  If  the  proportionality 
constant  is  nonnegative,  increase  z  until  x'jj=0  for  some  (i,j)  in  the  cycle;  there  will  be 

some  such  (i,j)  because  the  numbers  6jj  are  proportional  to  z  with  nonzero  proportionality 
constants  that  alternate  in  sign.  Otherwise,  decrease  z  to  achieve  the  same  end  Since  f()  is 
decreasing  in  each  variable,  the  x'  solution  is  at  least  as  good  as  the  x  solution,  with  one 
less  positive  variable.  By  repeating  this  operation  if  necessary,  an  optimal  solution  with  no 
cycles  can  be  found.  QED 

Since  an  optimal  conservative  solution  with  no  cycles  exists,  it  makes  sense  to  search 
for  the  optimal  solution  amongst  solutions  with  that  property.  A  major  advantage  of  such 
solutions  is  that  there  can  be  at  most  n+m-1  positive  variables  xjj.  To  see  this,  consider  the 
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nonoriented  bipartite  graph  G  formed  by  connecting  m  "row  nodes"  to  n  "column  nodes", 
with  i  being  connected  to  j  if  and  only  if  xjj>0.  Such  a  graph  must  consist  of  a  number  of 
connected  components  Ti,. .  .Tk,  each  of  which  is  not  connected  to  any  of  the  others. 
Some  of  these  components  may  contain  a  single  node  and  no  edges,  but  every  node  is  in 
exactly  one  component.  Since  each  component  is  connected  and  has  no  cycles,  it  is  a  tree 
(it  is  convenient  to  refer  to  even  singleton  components  as  "trees"),  so  G  is  a  "forest".  If  T\ 
has  rjc  nodes,  it  must  have  exactly  r^-l  edges  (Berge(1962)),  even  if  rk=l.  Therefore  the 
total  number  of  edges  is  n+m-K,  which  is  at  most  n+m-1. 

The  graph  G  described  above  will  be  called  the  "basis"  of  any  conservative  solution 
with  no  cycles.  More  generally, 
Definition:  Any  nonoriented  bipartite  graph  G  with  the  following  properties  is  a  basis: 

1)  The  nodes  of  G  are  all  of  the  m+n  row  and  column  nodes. 

2)  All  row  nodes  I  for  which  bi>0  have  at  least  one  incident  edge. 

3)  If  E(i,j)=0,  then  there  is  no  edge  (i,j). 

4)  There  are  no  cycles. 

The  requirement  that  a  solution  should  have  no  cycles  is  sufficient  to  establish  that  the 
corresponding  basis  is  a  forest,  but  it  is  not  sufficient  to  establish  a  one-to-one  relationship 
between  bases  and  solutions.  In  order  to  establish  such  a  relationship,  further 
assumptions  about  f()  are  needed. 
Definition:  f()  is  good  if  it  has  the  following  properties: 

1)  f()  is  differentiable  and  decreasing  on  the  interior  of  the  convex  set  S  described 
earlier.  Let  f  ()  be  the  gradient  of  f(). 

2)  There  exists  a  unique  function  g((i)  taking  values  in  S  such  that  f  (g((i))+(i=0  for 
all  |J>0,  where  by  (!>0  we  mean  |!j>0  for  j=l,...,n. 

3)  If  C  is  a  nonempty  subset  of  { l,...,n},  if  |i>0,  and  if  scalar  p>0,  then  the  equation 
^|ijgj(zjj.)=p  has  a  unique  positive  solution  z.  Furthermore  z  is  a  continuous  function  of 

j£C 

p  for  p>0. 

In  practice  the  methods  introduced  below  will  be  attractive  only  if  the  solution  z 
required  in  the  last  part  is  easily  computed,  since  the  equation  will  have  to  be  solved  many 
times.  Here  are  two  examples  of  good  functions: 

n 

Example  1:  f(y)=^  Vjexp(-yj),  with  S=Rn.  This  is  the  Search  Theory  function 
introduced  earlier,  except  that  pj  has  been  incorporated  into  Vj.  In  this  case  gj(|i)=ln(Vj/p:j) 

I>jgj(L0-p¥2>j  • 

JeC  J  LjeC     „ 


for  j=l,...,n,  and  ln(z)= 

Example  2:  f(y)=£  l/(l+yj),  with  S={y:  yj>-l  for  j=l n}.  Here  gJ(P-)=(Vv^J  -l) 


and  Vz= 


such  that 


I 

jeC 


j=l 


VM 


(p+Xki 


,  well  defined  for  p>0. 


The  Kuhn-Tucker  conditions  for  PI  are  that  there  should  exist  ?i,|i,x,  and  y 


5/5/89 


KT1)    Xi>jijE(i,j)  for  all  i,j 
KT2)    xjj>0  for  all  i,j 

KT3)    Xi=|ijE(ij)  when  xjj*0 

m 

KT4)    XxuE(i0)=yjforallJ 

i=l 
n 

KT5)    X  xij=bj  for  all  i,  and 
KT6)    y=g(n). 

The  object  now  is  to  use  the  equality  parts  of  the  KT  conditions  (KT3-KT6)  to  set  up  a 
one-to-one  correspondence  between  bases  and  solutions  for  good  objective  functions,  after 

which  solutions  corresponding  to  bases  will  be  called  "basic."  The  vectors  X  and  \i  will 
each  be  called  "multipliers."  Basic  solutions  will  have  xjj=0  for  nonbasic  edges  (i,j),  but 
the  following  theorem  permits  one  exception  to  support  the  pivoting  ideas  of  section  4. 

Theorem  2:  For  any  tree  T  let  R(T)  and  C(T)  be  the  row  nodes  and  column  nodes  of  T, 

respectively.  Let  G  be  a  basis  composed  of  trees  (Ti,. .  .,Tk),  and  suppose  that  ?ii=|!jE(i,j) 
for  all  edges  (i,j)  of  G.  Then 

a)  If  (x,y)  solve  KT4-KT5,  and  if  all  nonbasic  allocations  except  for  xtj  are  zero,  then 
for  any  component  T  of  G, 


X  M-jYj    -WxijE(I,J)5(C(T),J)  =     £  Jtibi 


jeC(T) 


-Xixu5(R(T),I),  (1) 


_  ieR(T) 

where  now  5(W,w)  is  1  if  w  is  in  the  set  W,  otherwise  0.  Furthermore, 

b)  If  (k,\x,y)  is  such  that  (1)  holds  for  every  component  T  of  G,  then  KT4-KT5  have 
a  unique  solution  x  for  which  all  nonbasic  allocations  except  xrj  are  zero. 

Proof  of  part  a)  Multiply  both  sides  of  KT4  by  ^j  and  sum  to  obtain 

m 

X   X  MjXijE(ij)=  X  M-jYj.  (2) 

jeC(T)  i=l  jEC(T) 

Since  xjj=0  when  jeC(T)  unless  either  i=I  or  ieR(T), 

X     Z  M-j>^ijE(i,j)  +  WxijE(I,J)5(c(T),J)  =  X  W  (3) 

jeC(T)  ieR(T)  jeC(T) 

But  |XjE(i,j)=A,i  when  jeC(T)  and  ieR(T)  by  assumption,  so 

X  Xj     X  XU   +  |iJxuE(I,J)5(c(T),J)  =  X  lijYj-  (4) 

i£R(T)        jEC(T)  jeC(T) 

Since  KT5  holds, 

X  xij   =bi  -  xu5({i},I)  for  all  ieR(T).  (5) 

jeC(T) 
Part  a)  follows  upon  substituting  (5)  into  (4). 
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Proof  of  part  b):   To  avoid  unnecessary  complication  assume  xtj=0;  otherwise  bj  and  yj 
can  be  adjusted  by  subtracting  xy  and  xtjE(IJ),  respectively.  Consider  any  component  T 
of  G.  We  must  show  that  there  is  exactly  one  way  to  assign  allocations  xjj  to  the  edges  <>t 
T  that  is  consistent  with  KT4-KT5.  This  is  trivial  if  T  is  a  singleton;  otherwise,  since  T  is 
a  tree,  there  must  be  at  least  1  pendant  (connected  to  exactly  one  other  node)  node.  If  this 
pendant  node  is  row  P,  let  row  P  be  connected  to  column  Q.  Then  xpQ  is  the  only  nonzero 
allocation  in  row  P,  and  is  therefore  necessarily  bp  by  KT5.  Define  y'  by 
y'Q=yQ-  bpE(P,Q),  y 'j=yj  for  j?^Q,  and  let  T'  be  the  tree  remaining  after  node  P  and  edge 

(P,Q)  are  deleted.  Since  ?ip=|J.qE(P,Q),  |iQy'Q=|iQyQ-?ipbp.  Therefore 

X  ^jy'j  =    X  ^'b»-  H" tne  pendant  node  is  column  Q,  let  it  be  connected  to  row  P. 
jeC(T)  ieR(T) 

Since  xpq  is  the  only  positive  allocation  in  column  Q,  necessarily  xpq  is  yQ/E(P,Q) 
according  to  KT4.  xpq  is  well  defined  because  E(P,Q)>0.  Let  b  p=bp-xpQ,  otherwise  let 
b'i=bj  for  i^P,  and  let  T  be  the  tree  resulting  from  deleting  column  Q  and  edge  (P,Q). 

Since  Xp=(IqE(P,Q),  ^pb'p=?tpbp-|J.QyQ.  Therefore    X  Kjvj  =    X  ^ib'i.    Since  in 

jeC(T')  ieR(T') 

either  case  T'  has  one  less  node  than  T,  and  since  the  proposition  is  true  for  graphs  with 
one  node,  the  theorem  is  proved  by  induction.  QED 

Multipliers  for  which  ?ij=}ijE(i,j)  on  basic  edges  are  easy  to  determine  if  there  is  no 

requirement  for  equation  (1)  to  hold.  The  reason  is  that  if  ^i  or  |ij  is  determined  at  any 

node,  then  the  same  thing  is  true  of  all  neighbors  of  the  node.  Thus  X  and  \i  can  be 
determined  in  any  component  T  of  G  by  assigning  an  arbitrary  positive  value  to  any  node. 

Let  (A*,|J*)  be  such  multipliers.  Then  any  other  set  (X,\i)  with  the  same  property  is  some 

scalar  multiple  of  (?i*,|i*).  Let  z>0  be  the  multiple.  Then  z  must  be  determined  so  that  (1) 

holds.  After  substituting  gj(z|i*)  for  yj,  (1)  is: 

X  |i*jgj(z^*)  -li*jxijE(I,J)5(C(T),J)=    X  ^*ib'  -  X*pcn6(R(T)J).  (6) 

jeCCT)  i£R(T) 

As  long  as  xij<bi,  the  right  hand  side  of  (6)  is  nonnegative.  Equation  (6)  will  then  have  a 
unique  positive  solution  z  as  long  as  f()  is  good.  Given  z,  KT3-KT6  have  the  unique 

solution  (k,\i)=z(X*,[i*)t  y=g(|i),  and  x  the  unique  solution  of  KT4-KT5.  The  x  part  is 
the  desired  basic  solution  corresponding  to  the  basis  G-hereafter  the  G-basic  solution. 
Thus  the  unique  G-basic  solution  can  be  easily  determined  as  long  as  xjj<bi . 

A  G-basic  solution  need  not  satisfy  KT2,  but  if  it  does  then  it  will  be  called  a  basic 
feasible  solution  and  G  will  be  called  "feasible".  Some  basic  feasible  solution  is  optimal, 
so  PI  could  be  solved  in  principle  by  examining  all  bases.  In  analogy  with  the  Simplex 
Method,  however,  the  object  is  to  avoid  this  by  utilizing  a  procedure  that  is  directional  in 
the  sense  of  considering  a  sequence  of  feasible  bases,  each  of  which  is  better  than  the  last. 
Two  of  these  will  be  considered  below.  In  each  case  we  will  prove  only  that  the  objective 
function  sequence  is  nonincreasing,  rather  than  decreasing,  with  equality  being  possible  in 
degenerate  cases.  Rather  than  deal  with  the  distracting  issue  of  degeneracy  here,  we  simply 
assume  nondegeneracy  for  the  next  two  sections.  There  is  every  reason  to  expect  that  the 
same  techniques  that  handle  degeneracy  in  the  Simplex  Method  will  also  work  here. 
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3 .    The  Manifold  Suboptimization  Method 

In  this  section  nonbasic  variables  will  always  be  0. 
A  G-basic  solution  is  the  optimal  solution  of  P2(G): 

mimimize  f(y) 


m 


subject  to         £  E(i,j)xij=yj;      l<j<n, 


i=l 

n 


£  xjj<bi;  l<i<m, 

xy=0  unless  (ij)EG. 

To  prove  this  note  that  the  G-basic  solution  solves  KT3-KT6.  But  KT3-KT6  are  the 
Kuhn  Tucker  conditions  for  P2(G),  and  the  Kuhn  Tucker  conditions  are  sufficient  for  a 
global  minimum  because  f()  is  convex  as  a  function  of  x.  P2(G)  does  not  require 
allocations  to  be  nonnegative,  so  the  feasible  space  is  not  compact.  Absent  the  constraints 
forcing  nonbasic  variables  to  be  0,  P2(G)  might  lack  a  (finite)  optimal  solution. 
Nonetheless,  the  G-basic  solution  is  optimal. 

Given  a  basis  G  and  a  corresponding  solution  x  feasible  in  PI,  suppose  KT1  is  not 
satisfied  for  some  (I,J).  The  natural  computational  step  is  to  simply  put  (I, J)  into  G.  Any 
cycle  that  forms  (a  cycle  will  form  if  and  only  if  I  and  J  are  in  the  same  tree)  can  always  be 
profitably  destroyed  by  removing  some  other  edge  (Kj,Li)  determined  as  in  Theorem  1. 
(I, J)  would  be  (i i  ,j  i )  in  that  theorem,  and  z  would  be  positive  because  the  reduced  cost  for 
(I,J)  is  negative.  Let  G1  be  G  with  (I,J)  added  and,  if  necessary,  (Ki,Li)  deleted.  Let  x1+ 
be  x  if  (Ki,Li)  is  not  deleted,  or  otherwise  the  solution  x'  of  Theorem  1.  x1+  is  better  than 
x  and  feasible,  but  is  not  in  general  the  G-basic  solution.  The  latter  solution  (call  it  x1)  is 
better  than  x1+  because  it  is  optimal  in  P2(G1)  while  x1+  is  merely  feasible,  x1  is  therefore 
basic  and  better  than  x.  However,  x1  may  not  be  feasible  in  PI. 

If  x1  has  one  or  more  negative  components,  let  x2+=x(l-a)+x1oc,  where  cc£[0,l]  is 
selected  so  that  x2+  is  the  closest  nonnegative  point  to  x1.  Let  (K2,L2)  be  an  edge  in  G1 
whose  allocation  is  driven  to  0  at  x2+.  Since  f()  is  convex  and  x1  is  better  than  x,  x2+  is 
also  better  than  x.  x2+  is  also  feasible  in  P2(G2),  where  G2  is  the  basis  obtained  by 
deleting  edge  (K2J-2)  from  G1.  Let  x2  be  optimal  in  G2,  so  that  x2  is  better  than  x2+.  If 
x2  also  has  negative  components,  continue  to  define  x3+,G^,  x^,  x4+,  etc.,  until  finally  xk 
is  nonnegative,  as  must  eventually  happen  because  an  edge  is  deleted  from  the  basis  at  each 
step.  xk  is  then  basic,  feasible  (in  PI),  and  (barring  degeneracy)  better  than  x.  Since  there 
are  only  finitely  many  bases,  none  of  which  can  repeat  because  the  objective  function 
decreases  at  each  step,  the  optimum  solution  must  eventually  be  found  in  finitely  many 
steps. 

The  above  algorithm  is  sufficiently  similar  to  ZangwiU's  (1970)  method  of  Manifold 
suboptimization  that  the  name  has  been  retained,  but  there  is  a  difference.  The  difference  is 
that  in  ZangwiU's  method  a  G-basic,  feasible  solution  that  is  not  optimal  would  lead  to  next 
considering  a  modified  basis  G'=G+(I,J),  where  (I, J)  is  some  edge  at  which  the  reduced 
cost  is  negative.  Unfortunately  G'  may  have  a  cycle,  in  which  case  P2(G')  may  not  have  a 
finite  optimal  solution.  The  method  described  above  works  because  every  manifold 
considered  has  a  basis  that  is  a  forest  of  trees. 

Manifold  suboptimization  essentially  operates  by  increasing  xtj  so  much  that  the  next 
basic  solution  may  not  be  feasible,  and  then  fixing  the  problem.  If  only  feasible  solutions 
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are  to  be  considered,  then  a  test  for  the  amount  that  xtj  can  be  increased  without  causing  an 
infeasibility  must  be  found.  This  is  the  subject  of  the  next  section. 

4 .    The  Extended  Simplex  Method 

In  this  section  it  will  be  necessary  to  consider  situations  where  nonbasic  variables  are 
temporarily  positive,  as  in  the  Simplex  Method.  Assume  that  G  is  a  feasible  basis,  but  that 
the  corresponding  basic  feasible  solution  is  not  optimal  because  KT1  is  not  satisfied  for 

some  (I,J).  Then  the  reduced  cost  Xt-|!jE(I,J)  is  negative,  so  xtj  should  be  increased  as  in 
Manifold  suboptimization.  Since  the  generality  will  be  needed  later,  suppose  that  the 
nonbasic  variable  xtj  is  fixed  at  some  level  b  between  0  and  bj. ,  and  that  the  balance 
equation  (6)  is  satisfied  for  each  component  T  of  G  by  z=l  when  xjj=b.  If  the  nonbasic 

variable  xjj  is  increased  from  b  to  b+A,  let  the  solution  of  (6)  for  z  be  Z(T,A).  If  f()  is  a 

good  function,  Z(T,A)  is  a  continuous,  positive  function  of  A  with  Z(T,0)=1.  The 

corresponding  multipliers  (^(A),|i(A))  are  the  same  as  (?i*,|i*)  except  that  multipliers  for 

nodes  in  T  are  multiplied  by  Z(T,A).  Let  y(A)=g(|i(A)),  let  x(A)  be  the  unique  solution  of 

KT4-KT5  where  only  basic  variables  are  nonzero  except  that  nonbasic  xjj=b+A,  let  A'  be 

the  largest  A  such  that  x(A)>0,  and  let  xkl(A')=0.  In  other  words,  basic  variable  (K,L)  is 

the  one  first  driven  to  0  by  increasing  A.  If  x(A)>0  for  all  A>0,  let  A'  and  (K,L)  be 
undefined.  If  (K,L)  is  defined  and  if  its  deletion  from  the  basis  breaks  any  cycle  that 
introducing  (I,J)  might  form,  then  the  next  basis  simply  replaces  (K,L)  with  (I,J).  (K,L)  is 
not  necessarily  defined  or  part  of  a  cycle  even  if  it  is  defined,  but  Algorithm  A  below 
accounts  for  these  possibilities.  The  input  to  the  algorithm  is  a  feasible  basis  G°, 
associated  multipliers  (.\°,ji°),  and  a  nonbasic  edge  (I,J)  with  xjj=0  and  negative  reduced 
cost.  The  output  is  a  different  feasible  basis  G'  at  least  as  good  as  G°,  together  with 

associated  multipliers  Ck',]!').  The  basic  idea  of  the  algorithm  is  to  increase  xtj  in  steps 
until  finally  the  reduced  cost  is  zero.  If  I  and  J  are  in  the  same  tree,  xtj  is  increased  in  step 
A4  until  some  edge  is  deleted.  Eventually  this  operation  will  put  I  and  J  into  different  trees. 

When  I  and  J  are  in  different  trees,  either  an  edge  is  deleted  from  one  of  the  trees  (A'<A  in 
step  A5)  or  else  (I, J)  is  finally  added  to  the  basis  in  the  terminal  step  A9.  Here  is 
Algorithm  A. 

Al)Letk=0,  XM) 

A2)  Let  \*=\k,  (i*=^ik,  G=Gk,  and  b=Xk. 

A3)  Let  Ti  be  the  tree  in  Gk  that  includes  I  and  Tj  be  the  tree  that  includes  J.  If  Tt>Tj 
go  to  A5. 

A4)  Let  T=Tj.  Determine  A'  by  solving  (6)  for  Z(T,A)  and  proceeding  as  in  the 
paragraph  above.  Let  (K,L)  be  the  basic  edge  in  T  whose  variable  is  first  driven  to  0,  let 
Gk+l=Gk-  (K,L),  (?lk+1,^lk+1)=(?l(A'),^!(A•)),  and  Xk+1=Xk+A'.  Go  to  A2. 

A5)  Let  T=Ti.  Determine  A'  by  solving  (6)  for  Z(T,A)  and  proceeding  as  in  the 
paragraph  above.  Let  (K,L)  be  the  basic  edge  in  T  whose  variable  is  first  driven  to  0.  Let 

(Ki,Lt>(K,L)  and  Aj=A'. 

A6)  Let  T=Tj.  Determine  A'  by  solving  (6)  for  Z(T,A)  and  proceeding  as  in  the 

paragraph  above.  If  (K,L)  is  undefined  let  Aj=°°;  otherwise,  let  (Kj,Lj)=(K,L)  and 
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Aj=A'. 

A7)  Let  A*  be  the  smallest  nonnegative  solution  of  the  equation 
?i*iZ(Ti,A)=|!*jZ(Tj,A)E(I,J),  with  A*=°<=  if  there  is  no  such  solution. 

A8)  Let  A=min{Ai,Aj,A*}.  Let  ^ik+1=^ikZ(Ti,A*)  for  ieTL  or?iik+i=XikZ(Tj,A*) 
for  ieTj,  or  otherwise  kjk+1=?ijk.  Define  )ik+1  similarly.    Let  Xk+1=Xk+A'. 

A9)  If  A'=A*,  set  G'=Gk+(U)  and  (?i',n')=(?ik+1,nk+1).   Stop. 

A 10)  If  A'=Ai,  let  Gk+1=Gk  -  (Ki,Li)  and  k=k+l.  Go  to  A2. 

All)  If  A'=Aj,  let  Gk+1=Gk-  (Kj,Lj)  and  k=k+l.  Go  to  A2. 

Theorem  3:  Algorithm  A  stops  after  finitely  many  steps.  The  output  is  a  feasible  basis 
different  from  the  input,  with  associated  objective  function  at  least  as  good. 

Proof:  A'  and  (K,L)  are  well  defined  in  A4  and  A5  because  row  I  is  included  in  Ti  and 

therefore  A'  cannot  exceed  bi.  If  I  and  J  are  included  in  the  same  subtree,  then  the  reduced 

cost  is  negative  when  A=0  and  proportional  to  Z(T,A),  therefore  negative  for  all  A.  If  I 

and  J  are  in  different  trees,  then  the  reduced  cost  is  negative  for  A<A*,  since  the  functions 

in  A7  are  continuous  in  A.  Since  A'<A*  in  A8,  the  reduced  cost  is  negative  throughout  and 
there  is  no  possibility  that  the  objective  function  might  increase.  It  is  also  clear  that  G'  is 
different  from  G,  since  G'  includes  (I,J)  and  G  does  not.  The  only  question  is  whether  A9 
is  encountered  after  finitely  many  steps. 

When  an  edge  is  deleted,  the  effect  on  the  basis  is  always  that  there  is  one  more  tree 
component  and  one  less  edge.  Since  A2  can  only  be  revisited  after  deleting  an  edge,  the 
number  of  edges  in  G°  is  an  upper  bound  on  the  number  of  times  A2  is  encountered.  Since 
A2  is  encountered  as  frequently  as  any  other  step,  the  proof  is  complete.       QED 

The  pivoting  operation  described  above  differs  from  the  comparable  Linear 
Programming  operation  in  that  the  number  of  edges  that  have  to  leave  the  basis  in  order  to 
get  (I,J)  in  is  not  always  1;  it  may  be  any  nonnegative  integer.  It  differs  essentially  in  this 
respect  from  Zangwill's  (1970)  Convex  Simplex  Method  (CSM).  In  the  CSM  the  number 
of  basic  variables  is  constant,  and  every  pivot  replaces  a  basic  variable  with  a  nonbasic 
variable.  When  a  local  minimum  is  encountered  in  CSM,  the  corresponding  variable 
remains  nonbasic,  but  positive,  and  consequently  there  is  no  unique  solution  associated 
with  a  given  basis.  In  the  method  discussed  above,  the  number  of  basic  variables  is  not 
constant,  but  the  idea  that  nonbasic  variables  are  zero  is  preserved,  along  with  the  one-to- 
one  relationship  between  bases  and  solutions.  This  relationship  implies  that  the  number  of 
pivots  is  bounded  by  the  number  of  bases,  whereas  the  number  of  pivots  required  by  the 
CSM  may  be  unbounded. 

5 .    Implementation  in  the  Exponential  Case 

Two  FORTRAN  programs  have  been  written  that  implement  the  methods  of  sections  3 
and  4  in  the  Exponential  case:  MANIFO  in  the  case  of  Manifold  suboptimization  and 
SIMPLX  for  the  Extended  Simplex  Method.  Each  of  these  programs  exploits  the  forest 
structure  of  the  basis  to  minimize  storage  and  facilitate  computations.  The  forest  is  actually 
represented  as  a  single  tree  by  introducing  an  "earth"  node  to  which  all  of  the  components 
of  the  basis  are  connected  through  a  root  node  in  each  tree;  however,  in  this  exposition  the 
term  "tree"  will  continue  to  mean  a  component  of  the  basis  as  heretofore.  Most  of  the 
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operations  associated  with  pivoting  are  similar  to  those  encountered  in  transshipment 
problems.  Since  these  are  well  discussed  in  Bradley,  Brown,  and  Graves  (1977),  only  the 
aspects  associated  with  the  nonlinear  nature  of  the  objective  function  will  be  discussed 
here. 

Consider  first  the  problem  of  deleting  an  edge  from  the  basis.  If  the  edge  is  in  tree  T, 
then  deleting  it  amounts  to  pruning  off  a  branch  from  T.  A  reduced  version  of  T  will 
remain,  along  with  a  new  tree  T'  consisting  of  the  pruned  off  branch.  Although  equation 
(1)  will  hold  in  all  trees  before  the  pruning,  it  generally  will  not  hold  in  either  T  or  T' 

afterwards.  The  multipliers  (k,[i)  must  be  scaled  to  make  (1)  hold  in  T  and  T,  bearing  in 

mind  that  y  depends  on  |i  and  that  nodes  in  T  will  have  a  different  z-factor  than  nodes  in 
T.  These  z-factors  can  be  computed  using  the  formula  given  in  Example  1  of  Section  2. 

The  computation  in  T  requires  the  sum    2^  M-j-  This  number  could  be  computed  by 

jeC(T') 
summing  [ij  over  the  appropriate  set,  but  in  fact  both  MANIFO  and  SIMPLX  utilize  an 

additional  array  SDM(k)  to  store  the  sum  of  |ij  over  the  set  consisting  of  k  and  all  of  its 
column  successors  in  whatever  (rooted)  tree  k  belongs  to.  Thus,  if  (I,J)  is  deleted  from  a 
tree  with  root  node  K,  and  if  (say)  J  is  the  predecessor  of  I,  then  the  required  sum  in  T  is 
SDM(I)  and  the  required  sum  in  T  is  SDM(K)-SDM(I).  Calculating  the  z-factor  then 

requires  a  single  exponentiation  in  each  tree.  After  updating  (k,\i,y),  the  allocations  in  T 
and  T'  are  computed  as  described  in  part  b)  of  Theorem  2.;  a  preorder  traversal  array 
permits  this  to  be  done  in  one  pass.  The  allocations  themselves  are  stored  in  array  X,  with 
X(k)  representing  the  allocation  from  node  k  to  its  predecessor  if  k  is  a  row,  or  from  the 
predecessor  to  k  if  k  is  a  column.  All  arrays  are  of  length  n+m+1  in  both  programs  (node 
length);  in  fact,  the  only  function  requiring  mn  operations  is  that  of  determining  the  entering 

edge  by  searching  for  the  largest  ratio  |ijE(i,j)Ai . 

The  addition  of  edge  (I,J)  to  the  basis  requires  the  joining  of  one  tree  to  another, 
possibly  after  or  in  conjunction  with  another  pruning  operation.  If  T  containing  J  is  to  be 
joined  to  T,  T'  is  first  "rehung"  so  that  the  root  node  is  J.  Next  T'  is  scaled  so  that 

?q=|ijE(I,J),  which  requires  a  logarithm  in  updating  the  potentials  in  T.  The  unified  tree  is 
then  scaled  as  in  the  previous  paragraph.  It  is  during  the  subsequent  calculations  of  X  that 
the  possibility  of  negative  allocations  must  be  provided  for  in  MANIFO. 

SIMPLX  uses  all  of  the  arrays  in  MANIFO  plus  one  additional.  The  additional  array 

W  is  associated  with  determining  how  big  A  can  be  before  some  basic  variable  becomes  0, 

as  will  be  explained  shortly.  In  the  Exponential  case  the  function  Z(T,A)  satisfies 

ln(Z(T,A))=D(T)A,  where 

D(T)=[^*5(T,I)-|ij*E(I,J)6(TJ)]/  £  M-j-  <7) 

JEC(T) 

Note  that  D(T)  is  easy  to  compute  because  the  denominator  is  already  stored,  and  also  that 
D(T)  will  be  negative  in  step  A5,  positive  in  A6,  and  negative  (since  the  reduced  cost  is 

negative)  in  A4.  Multiplying  (ij  by  Z(T,A)  is  equivalent  to  subtracting  ln(Z(T,A))  from  yj, 

so  yj(A)=yj(0)-D(T)A  for  jeC(T).  Since  this  is  a  linear  function  of  A,  allocation  xjj(A)  will 

also  be  a  linear  function  of  A  for  all  edges  in  T.  The  array  element  W(k)  contains  the  rate  at 

which  the  allocation  associated  with  node  k  decreases  with  A.  The  calculation  of  W  is 
equivalent  in  difficulty  to  the  calculation  of  X  in  MANIFO.  Given  W,  it  is  simple  to  obtain 
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the  limiting  value  A'  by  comparing  X/W  ratios.  The  X  array  can  be  updated  once  A'  is 
known. 

Both  procedures  require  a  basis  for  a  starting  point.  This  is  obtained  by  first  including 
all  edges  of  the  form  (k,k).  If  there  are  more  columns  than  rows,  this  will  leave  some  trees 
consisting  of  singleton  columns.  If  there  are  more  rows  than  columns,  then  row  n+k  is 
assigned  to  column  k  for  l<k<n,  etc.,  until  finally  each  row  has  been  included  in  some 
edge.  The  resulting  basis  will  be  feasible,  but  unfortunately  some  edge  (i,j)  for  which 
E(i,j)  =0  may  be  included.  The  simplest  (and  the  operative)  remedy  is  to  simply  change 
E(i,j)  to  some  small  positive  number,  since  in  any  case  either  procedure  will  quickly  delete 
such  an  edge. 

There  is  potential  for  roundoff  errors  to  accumulate  in  the  repeated  updating  of  floating 
point  arrays,  so  MANIFO  and  SIMPLX  both  use  double  precision  arithmetic.  Neither 
program  incorporates  any  protection  against  degeneracy,  but  so  far  there  have  been  no 
instances  of  cycling. 

To  give  the  reader  a  better  idea  of  the  actual  operation  of  the  algorithms,  Figure  1 
shows  the  sequence  of  bases  considered  by  MANIFO  in  the  process  of  solving  a  problem 
where  m=3  and  n=4.  The  basis  is  shown  as  a  single  tree  with  row  nodes  being 
distinguished  from  column  nodes  (as  they  are  in  MANIFO)  by  representing  row  nodes 
with  negative  numbers.  Node  0  is  "earth".  With  the  convention  of  including  all  edges 
incident  to  node  0,  a  basis  always  has  exactly  n+m  edges.  Although  X  is  stored  as  a  node 
length  array  in  MANIFO,  the  allocations  are  shown  in  figure  1  in  the  more  familiar  matrix 
form.  There  are  (24— 1)3=3375  bases  for  this  problem,  of  which  MANIFO  considers  8 
feasible  and  2  nonfeasible  before  finding  one  that  is  optimal.  In  the  course  of  doing  this 
MANIFO  requires  18  logarithms  or  exponentials.  These  statistics  are  the  same  for 
SIMPLX  in  this  small  problem;  differences  emerge  only  in  problems  where  MANIFO 
introduces  more  than  1  negative  allocation  at  a  time. 

6 .    Computational   Comparisons 

MANIFO  and  SIMPLX  were  compared  on  randomly  generated  problems  to  determine 
which  program  is  fastest.  Problems  generated  were  such  that 
E(ij)  is  uniform  [0,1] 
bj  is  uniform  [0,10] 
V(j)  is  uniform  [0,100] 
Eleven  (m,n)  pairs  were  tested:   (10,10),  (7,20),  (20,7),  (20,20),  (30,30),  40,40), 
(20,40),  (40,20),  (50,50),  (60,60),  and  (70,70).  Each  of  eleven  (E,b)  configurations  was 
combined  with  each  of  five  V  configurations  to  generate  a  problem.  For  each  of  these  55 
problems,  the  number  of  nonlinear  operations  (logarithms  and  exponentials)  NL  was 
recorded  along  with  the  CPU  time  T  in  seconds  required  to  set  up  the  initial  guess  and 
solve  the  problem,  m,  n,  NL,  and  T  were  then  transformed  by  taking  logarithms  and  a 
linear  regression  performed. 

There  turns  out  to  be  very  little  difference  between  MANIFO  and  SIMPLX.  The  fitted 
formulas  are  in  each  case  approximately 

T  =  (.344sec)(m/46)1-05(n/122)-8exp(m/46+n/122)    (R2=99%)       (8) 

NL=1.4mVn     (R2  =  98%)  (9) 

All  runs  were  made  on  the  Naval  Postgraduate  School's  IBM3033AP.  Since 
GAMS/MINOS  is  also  available  on  that  system,  in  one  instance  a  problem  with  m=20  and 
n=20  was  solved  using  GAMS/MINOS.  The  time  required  was  1.88  seconds  ("MINOS 
TIME"  in  the  output),  which  is  47  times  as  large  as  the  time  required  by  MANIFO  (.04 
seconds)  for  the  same  problem.  Further  comparisons  were  not  made  because  MINOS  is  a 
general  purpose  solver  that  does  not  exploit  the  network  nature  of  the  constraints.  A  better 
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comparison  would  be  with  the  GENOS  network  optimizer  (Mulvey  and  Zenios  (1987;;. 
Attempts  to  accomplish  this  are  in  progress,  but  so  far  the  lack  of  an  exponential  function  in 
GENOS  1 .0  has  proved  to  be  a  roadblock. 

The  reason  for  keeping  track  of  the  number  of  exponentials  and  logarithms  NL  in 
running  MANIFO/SIMPLX  is  that  it  was  initially  anticipated  that  such  nonlinear  operations 
would  require  relatively  large  amounts  of  time.  This  turns  out  not  to  be  the  case.  When 
(m,n)  =  (70,70),  formulas  (8)  and  (9)  predict  T  =  .28  seconds  and  NL  =  820.  Since  an 
IBM3033AP  requires  only  about  4|isec  for  a  double  precision  logarithm  or  exponential,  the 
amount  of  time  spent  in  nonlinear  operations  is  orly  about  (820)(4xl06)  =  .0033  seconds; 
MANIFO  and  SIMPLX  are  each  mainly  occupied  with  manipulating  and  scaling  the 
various  arrays  required  to  represent  the  current  basic  solution  as  a  tree,  rather  than  solving 
the  balance  equation.  This  is  encouraging  for  applications  where  the  objective  function, 
although  "good,"  might  require  a  numeric  rather  than  analytic  solution  of  the  balance 
equation. 
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Figure  1 
Solution  of  a  3x4  Problem  Using  MANIFO  (Continued) 
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Figure   1. 
Solution  of  a  3x4  Problem  Using  MANIFO 
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